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The  overall  method  contains  two  parts,  a  solution  algorithm 
and  a  mesh  selection  algorithm.  The  solution  algorithm  is  a 
finite  element-Galerkin  method  on  trapezoidal  space-time  element 
using  either  piecewise  linear  or  cubic  polynomial  approximations 
and  the  mesh  selection  algorithm  builds  upon  similar  work  for 
variable  knot  spline  interpolation. 


A  computer  code  implementing  these  algorithms  has  been 
written  and  applied  to  a  number  of  problems.  These  computations 
confirm  that  the  theoretical  error  estimates  are  attained  and 
demonstrate  the  utility  of  variable  mesh  methods  for  partial 
differential  equations.. 
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The  overall  method  contains  two  parts,  a  solution  algorithm  and  a  mesh 
selection  algorithm.  The  solution  algorithm  is  a  finite  element-Galerkin 
method  on  trapezoidal  space-time  elements,  using  either  piecewise  linear  or 
cubic  polynomial  approximations  and  the  mesh  selection  algorithm  builds  upon 
similar  work  for  variable  knot  spline  interpolation. 

A  computer  code  implementing  these  algorithms  has  been  written  and  applied 
to  a  number  of  problems.  These  computations  confirm  that  the  theoretical  error 
estimates  are  attained  and  demonstrate  the  utility  of  variable  mesh  methods  for 
partial  differential  equations. 
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1.  Introduction 


In  this  paper  we  construct  an  adaptive  grid  finite  element  procedure  to  find 
numerical  solutions  of  M-dimensional  vector  systems  of  partial  differential  equa¬ 
tions  having  the  form 

(1.1)  Lu  :=  ufc  +  f(x,t,u,u^)  -  lD(x,t,u)u^]^  =  0,  0  <  x  <  a,  t  >  0, 
subject  to  the  initial  and  linear  separated  boundary  conditions 


(1.2) 


u(x,0)  =  u  (x)  ,  0  <  x  <  a, 


B  u(0,t] 


:=  A11(t)u(0,t) 


+  A, .  (t)u  (0,t)  =  b  (t) , 
-12  ~x  -1 


(1.3) 

B2u(a,t)  :=  A21(t)u(a,t)  +  A22 (t) u^  (a , t)  =  b2  (t) ,  t  >  0. 

There  are  initial  boundary  conditions  at  x  =  0  and  k2  terminal  boundary  con¬ 
ditions  at  x  =  a.  We  are  primarily  concerned  with  solving  diffusion  problems 
where  D  is  positive  definite  and  k^  +  k2  *  2M;  however,  we  will  not  restrict 
ourselves  to  this  case,  but  instead  we  assume  that  conditions  are  specified  so 
that  (1.1) -Cl. 31  has  an  isolated  solution. 

Problems  of  the  above  form  arise  in  many  applications  which  model  problems 
as  diverse  as  heat  conduction  (cf .  Friedman  {16j ) ,  determining  bacterial  motion 
(cf.  Keller  and  Odell  [25,30]),  combustion  (cf.  Kapila  [24]),  chemical  reactions 
(cf.  Fife  [14]),  population  dynamics  (cf.  Hoppensteadt  [20]),  and  convecting 
flows  (cf.  Batchelor  14] ) .  Therefore,  a  general  purpose  code  to  solve  (1.1)- 
(1.3)  numerically  would  be  extremely  useful. 

Many  of  the  problems  mentioned  above  have  solutions  which  contain  sharp 
transitions  such  as  boundary  layers,  shock  layers,  or  wave  fronts.  In  order  to 
resolve  such  nonuniformities  using  a  minimum  number  of  mesh  points  it  is  desir¬ 
able  to  concentrate  the  mesh  within  the  transition  layers.  Since  these  transi¬ 
tion  layers  can  move,  it  is  all  the  more  desirable  for  the  mesh  to  adapt  itself 
with  the  evolving  solution.  To  do  this  we  develop  methods  that  (i)  discretize 


(1.1)— (1-3)  on  a  nonuniform  mesh  and  (ii)  determine  the  proper  mesh  point  loca¬ 
tions. 


We  discretize  (1.1)- (1.3)  using  a  finite  element  Galerkin  method  on  trape¬ 
zoidal  space-time  elements.  This  approach  is  similar  to  that  of  Jamet  and 
Bonnerot  [6,22]  and  it  was  chosen  because  it  is  generally  easier  to  generate 
high  order  approximations  to  partial  differential  equations  on  a  nonuniform  mesh 
with  finite  element  methods  than  with  finite  difference  methods.  The  accuracy 
and  order  of  convergence  of  our  methods  are  analyzed  in  Davis  [11]  and  are 
demonstrated  experimentally  in  Section  4  of  this  paper. 

Adaptive  mesh  selection  strategies  typically  involve  some  recomputation  of 
the  solution.  That  is,  an  initial  solution  is  computed  on  a  coarse  mesh  and 
this  is  used  to  determine  whether  to  add  mesh  points  to  some  portion  of  the 
domain  and  redo  the  calculation,  redo  the  calculation  using  a  more  accurate 
method,  redo  the  computation  using  some  combination  of  these  methods,  or  accept 
the  present  computation.  Algorithms  of  this  general  type  have  been  developed 
and  successfully  applied  to  adaptive  quadrature  (cf.  eg.  Rice  [33]  and  Lyness 
and  Kaganove  [28]),  two-point  boundary  value  problems  (cf.  eg.  Childs  et  al 
[9]),  elliptic  boundary  value  problems  (cf.  eg.  Carey  [8]  and  Brandt  [7]),  and 
parabolic  and  hyperbolic  problems  (cf.  eg.  Berger  et  al  [5]  and  white  [37]). 

Primarily  because  of  the  expense  involved  in  recomputing  the  solution  of 
the  partial  differential  equations  at  possibly  every  time  step  we  have  developed 
an  algorithm  which  initially  places  a  fixed  number  of  mesh  points  in  optimal 
locations  and  then  attempts  to  move  them  so  that  their  locations  remain  optimal. 
Algorithms  of  this  type  have  been  used  by  Lawson  [26],  deBoor  [12,13],  and  Jupp 
[23]  for  variable  knot  spline  interpolation  and  it  is  their  work  that  motivated 
our  mesh  selection  algorithms. 

A  different  approach  to  this  problem  was  proposed  by  Miller  and  Miller  [29] 
and  later  extended  by  Galinas  et  al  [17] .  They  approximated  the  solution  of 
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parabolic  partial  differential  equations  by  piecewise  linear  polynomials  where 
both  the  polynomial  coefficients  and  the  mesh  on  which  they  were  defined  were 
unknown  functions  of  time.  These  functions  were  determined  by  minimizing  the 
least  squares  residuals.  They  found  that  the  mesh  points  could  coalesce  in 
certain  situations  and  they  avoided  this  by  adding  a  number  of  spring  and 
damping  terms  as  constraints  to  the  equations. 

One  advantage  of  the  above  approach  is  that  it  readily  extends  to  higher 
dimensional  problems.  However,  we  are  not  convinced  that  it  is  necessary  to 
couple  the  solution  and  mesh  selection  methods.  This  can  dramatically  increase 
the  size  of  the  discrete  system  without  offering  any  corresponding  increase  in 
order  of  accuracy.  Furthermore,  the  entire  solution  procedure  must  halt  if  an 
acceptable  mesh  cannot  be  calculated.  Under  the  same  circumstances  our  methods 
can  continue  to  compute  a  solution  on  a  suboptimal  mesh.  Since  both  methods  are 
under  development  we  have  not  attempted  any  detailed  comparisons. 

In  Section  2  of  this  paper  we  develop  a  finite  element  Galerkin  approxima¬ 
tion  to  (1.1) -(1.3)  using  trapezoidal  space-time  elements.  In  Section  3  we 
describe  a  practical  and  efficient  mesh  selection  procedure  that  approximately 
minimizes  the  L2  error  of  the  computed  solution.  In  Section  4  we  apply  the 
method  to  a  number  of  problems  and  discuss  the  computed  results.  Finally,  in 
Section  5  we  present  an  overall  discussion  of  this  effort  and  some  suggestions 


for  future  work. 
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2.  Finite  Element  Formulation 

We  discretize  (1.1,2, 3)  using  a  finite  element-Galerkin  procedure.  To  this 
end,  let  be  the  strip 

(2.1)  S  !=  {(x,t)|o  <  x  <  a,  t  <t<t  ,1, 

choose  "test"  or  "weight"  functions  v(x,t)EC° {S^) ,  multiply  (1.1)  by  v,  integrate 

over  S^,  and  integrate  the  time  derivative  and  diffusive  terms  by  parts  to  get 

tn+l  a 

F(u,v)  /  /  {-u  v  +  f (x , t , u ,u  ) v  +  D(x,t,u)u  v  }dxdt 

t0~t-  ~  ~x  ~  ~~xx 

(2.2)  n  t  t 

a  n+1  n+1  a 

+  /  u  v  dx|  -I  D  u  v  dt |  =  0. 

0  ~  t=t  t  ~X  x=0 

n  n 

Equation  (2.2)  is  called  the  Galerkin  form  of  the  problem  and  any  function  u  that 
satisfies  (2.1)  and  the  initial  and  boundary  conditions  (1.2,3)  is  called  a  "weak 
solution. " 

We  introduce  a  mesh  (O  =  x?  <  x"  <  ...  <  x"  =  a)  at  t  =  t  and  a  different  mesh 

12  N  n 

{O  =  <  x2+^  =  a)  at  t  =  ^n+i*  We  connect  the  corresponding  points  x^ 

and  x?+^  by  straight  lines  and  consequently  divide  the  strip  into  a  set  of 

N  -  1  trapezoids.  We  let  x.  (t)  denote  the  straight  line  connecting  xn  and  xn+1 

l  ii 

and  T°  denote  the  trapezoid  with  vertices  (x?,t  ),  (x?  , ,t  ) ,  (xn+^,t  ,), 

l  in  l+l  n  l+l  n+1 

<x"+1,tn+1)  (cf.  Figure  1). 

We  approximate  u(x,t)  on  S  by  U(x,t)€  U(S  )  which  has  the  form 
~  n  ~  K  n 

K 

(2.3)  U (x,t)  =  1  2i (t)^i (x't) • 

i=l 

The  "trial"  functions  tj^(x,t),  i  =  1,2,...,K,  can  be  used  to  construct  a  basis 

for  U„(S  ).  They  are  selected  to  be  of  class  C°(S  )  and,  in  finite  element 

methods,  to  have  small  support.  Particular  choices  of  4k  are  given  in  Section 

2.1;  herein,  it  suffices  to  note  that  4,  is  nonzero  only  on  T?  ,UTn  and  that 

i  1  l-l  i 

K  must  be  at  least  N. 
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We  determine  U  on  by  solving  a  discrete  problem  of  the  form 

(2.4)  P  U(x,0)  =  P  u°(x)  ,  n  =  0 

(2.5)  F  (U,v)  =  0,  Vv  el 

—  —  K 

<J-6>  ii  ■  ii'W-  ?2  "‘■■'W  ■  i/w- 

Here  V  is  a  finite  dimensional  space  of  C° (S  )  functions  that  depends  on  the 
K  n 

boundary  conditions  (cf.  Section  2.3),  P  is  an  interpolation  operator  (cf.  Section 

A  A  A  A 

2.3),  ,  b^,  B2»  b2  are  approximations  of  B^ ,  b^ ,  B2>  b2  obtained  by  numerical 

A 

integration  (cf.  Section  2.3),  and  F(U,v)  is  an  approximation  of  F(U,v)  obtained  by 

evaluating  the  integrals  in  (2.2)  numerically  (cf.  Section  2.2).  Equations  (2.5,6) 

result  in  an  MK  dimensional  nonlinear  algebraic  system  for  determining  the  Galerkin 

coordinates  c.  (t  ) ,  i  =  1,2, .. .K,  in  terms  of  c.  (t  ) ,  i  =  1,2, ...  ,K.  Since 
~i  n+i  ~ l  n 

c±  (0) ,  i  =  1,2, ... ,K,  are  determined  from  the  initial  conditions  (2.4),  equations 

(2.4-6)  define  a  marching  algorithm  for  determining  U(x,t)  in  successive  strips 

S  ,  n  =  0,1, . . .  . 
n 

If  there  were  no  boundary  conditions  we  would  select  (x,t),  i  =  1,2,...,K, 
as  a  basis  for  K  .  This  prescription  has  to  be  modified  slightly  for  i  =  1  and/or 

J\ 

i  =  K  (cf.  Section  2.3)  since  boundary  conditions  are  generally  imposed;  however, 
it  is  still  appropriate  to  write  F(U,v)  as  a  sum  of  contributions  from  each 
trapezoid.  Thus, 

N-l 

F (U,v)  =  E  If  (-U  v  +  f(x,t,U,U  ) v  +  D(x,t,U)U  v  dx  dt 
~  -  Tn  ~  t  -  ~x  -  ~  -x  x 

(2.7)  mi  v  (t)  t  t  a 

+  I  [/  1+1  U  v  dt]  n+1  -  [/  n+1  DU  v  dt]  =  O.VveV 

i“l  Xi(t)  "  t=tn  "n  “X  X=° 

Since  the  bases  for  both  the  trial  and  test  spaces  have  small  support  most 

of  the  integrals  in  (2.7)  will  be  zero.  The  algebraic  system  (2.5,6)  will  be 
sparse  and,  hence,  it  may  be  solved  effeciently. 
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2.1.  Selection  of  a  Basis 


A  simple  way  to  construct  a  basis  on  trapezoidal  elements  that  satisfies 

the  necessary  continuity  requirements  and  has  small  support  is  to  apply  a  local 

transformation  that  maps  each  trapezoid  onto  a  rectangle.  The  inverse  of  this 

transformation  on  T?  is 
x 

n  ,  ,  n  n.  .£+1.  ,  .  n+1  nv 

x  =  xi +  (xi+rxi)  <— 5 1  +  (*i  -xi)T 


(2.8) 


,  ,  n+1  n+1  n  n.  ,S+1. 

+  <Xi+1  -  X.  -  X.+1  +  X.)  (~ )T, 


t  +  (t  -t  ) T  . 
n  n+1  n 


It  maps  the  rectangle 


(2.9)  R  =  {(S,T)  |-  1  £  S  <_  1,  0  £  T  <  l) 

in  the  (S,T)  plane  onto  t"  in  the  (x,y)  plane. 

We  choose  this  basis  so  that  <Jk  (x,t)  is  a  function  of  S  only  on  .  To 
be  specific,  we  currently  allow  <t>^(x,t)  to  be  either  a  piecewise  linear  or  a 
piecewise  Hermite  cubic  polynomial  in  S  on  T^. 

For  piecewise  linear  approximations  we  construct  a  basis  in  terms  of  the 
canonical  basis  function 

A 

(2.10)  (MS)  =  (1-0/2 ,  -1  <  s  <  1, 


by  defining 
(2.11)  <(>i(x,t) 


<MS)  ,  (x,t)fi" 

A  X 

'  <M-S)  ,  (x , t)  e  T"_r  i  =  1,2 . K  =  N. 

0  ,  otherwise 


Thus,  the  dimension  of  the  trial  space  U  is  K  =  N.  Along  the  line  x.(T)  joining 

K  3 

n  ,  n+1  . 

x.  and  x.  we  have 
3  3 

(2.12)  $ . (x , (T ) ,  t (T ) )  =  6  .,  0  <  T  <  1, 


where  6  is  the  Kronecker  delta.  Using  (2.3)  this  implies  that 


(2.13) 


cilt)  =  Ui(T)  :=  U(x  (T)  ,  t  (T  ) ) 


Thus,  since  only  and  +  1  are  nonzero  on  T?  we  have 


(2.14) 


U(x,t)  *U.(T)<M£)  +  U.  ,  (T)^(-C),  (x,t)CT. 
~  -1  -1+1  1 


For  piecewise  cubic  Hermite  approximations  we  construct  a  basis  in  terms  of 
the  two  canonical  basis  functions 

(2.15)  <MS)  =  j(i-02<2+£) .  \(K)  «  j(i+0  (K)2,  -1  <  K  <  1, 


by  defining 

*  /v 

Ha  ,  (x,t)  c  t? 

(2.16a)  <i’2i-l(X,t)  =<  (*,t)  (Tj  ,  i  =  1,2 . N 

0  ,  otherwise 


X(H  ,  (x,t)  c  tV 

(2.16b)  *2i(x,t)  (x,t)  C  t"_x,  i  =  1,2, ...  ,N 

0  ,  otherwise 


Thus,  the  dimension  K  of  the  trial  space  in  2N. 

We  note  that  <f_.  ,(x,t)(c\s  )  with 
2i~l  n 


(2. 17a,b)  <J>2i_1(x.  (T)  , t  (T ) )  *  6  ,  <t2i_1  (X.  (T)  , t  (T )  )  =  0,  0  <  T  <  1, 

1  3  13  x  ^ 

but  (J>  (x,t)  f  C°(S  )  with 

2i  n 

<2.17c,d)  <f>2i(x  (T),t(T))  =  0,  <f>2_  (x .  (T )  ,  t  (T ) )  =  6.  ,/xr(T)t  0  <  I  <  1. 
J  x  3  X-D  ^ 

The  function  x.  (T )  is  easily  computed  from  (2.8)  as 

,  ,  1 ,  n  n,  1  ,  n+1  n+1  n  ,  n.  . ,  ,  ^ .  ,  n 

(2.18)  x.(T)  =  -(x  -x.)  +  -(x  -x.  -x .  .  +x . )  T  ,  if  (x ,  t )  £  T . . 

23  +  13  2  3+1  3  3+1  3  3 


Thus,  x^(T)  is  different  on  each  trapezoid  (unless  the  mesh  is  uniform  and 
rectangular)  and  $2^  (x,t)  jumps  as  x  crosses  x^(T).  However,  using  (2.3)  and 
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(2.17)  we  can  make  U(x,t)  of  class  C1  (S  )  by  selectir 


(2.19a) 


~2i-l(t)  =  ~i(T)  :=  U(x.(T)  ,t(x))  , 


(2.19b) 


C  (t)  =  x  (T)U  (I)  :=  X  (T)  u  (X.  (T)  ,t(T)  )  , 

'•21  t,  ~X^  f,  -XI 


Thus,  on  T .  we  have 
1 


(2.20)  U  ( x ,  t )  =  Ui(T)'|(C)  +  Ui+i  +  Uy  (T)x  (T)x(^) 

i  ^ 

A 

-U  (T)X  (T)X(-C). 
l+l  ^ 

2.2  Numerical  Integration 

Ignoring  the  boundary  conditions  for  the  moment,  we  choose  v  =  d  (x,t) 
according  to  either  (2.11)  or  (2.16)  and  use  (2.8)  to  transform  (2.7)  to 


(2.21)  F(U,<|>.)  =  l  I.(U,<{>.)  -  I  (U,<f>  .)  =  0,  j  =  1,2 . K, 

“  "  ]  i=l  ~X  ~  2  ~  1 

•where 


(2.22a) 


I, (U,v)  =  /  /  {-U  v  £  +  f  (x,t,U,Ur£  ) v 
~  0~1  -  i,  t  ~  ~  -i,  x 


(2.22b) 


+  D (x , t , U) U  vr?2} | j|d£di  +  f1  U  v  x  d^l1 

~  ~  ^  X  -1  ~  K  T  =  0 

.  x=a,£=l 

ID(U,v)  =  f  D(x,t,U)Urv£  t  dl| 

~B  '  0  ~  XT  x=0,£=-l 


The  functions  £^,  £^,  x^i  t^ ,  and  |j]  the  Jacobian  of  the  transformation 
can  be  computed  from  (2.8). 

In  order  to  complete  the  specification  of  our  numerical  method  we  need  to 
select  quadrature  rules  for  evaluating  the  integrals  in  (2.22).  We  use  the 
Trapezoidal  rule  to  evaluate  the  T  integrals  and  a  three  point  Gauss-Legendre 
rule  (cf.  Abromowitz  and  Stegen  [1],  Chap.  25)  to  evaluate  the  £  integrals. 

The  latter  was  chosen  because  it  is  known  (cf.  Strang  and  Fix  [35])  to  have  the 
same  order  discretization  error  as  our  finite  element  method  with  cubic 
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approximations  and  the  exact  integration  of  (2.22).  At  present,  we  also  use  the 
three  point  Gauss-Legendre  rule  for  linear  approximations  although  it  is  more 
accurate  than  necessary  in  this  case  and  therefore  somewhat  inefficient. 

Upon  use  of  the  above  mentioned  quadrature  rules  equations  (2.25)  become 
^  N-l  „ 

(2.23)  F  (U,<f> . )  =  Z  I.(U,<f>.)  -  ID(U,<1>.)  =  0,  j  =  1,2,... ,K, 

~  ~  :  i=i  ~1  ~  3  ~b  ~  j 

where  1.  and  1^  denote  the  approximations  of  (2.22)  that  are  obtained  by  numerical 
integration. 

2.3.  Initial  and  Boundary  Conditions.  Solution  Technique, 

The  solution  U  is  determined  on  by  solving  (2.23)  together  with  the 
initial  and  boundary  conditions  (2.4)  and  (2.6),  respectively.  We  satisfy  the 
initial  conditions  (and  implicitly  define  the  interpolation  operator  P  of  (2.4)) 
by  requiring 

(2.24a)  U°  =  u°  ( x i ) ,  i  =  1,2 . N, 

for  both  linear  and  cubic  approximations,  and  additionally 


(2.24b)  U°  =  u°(x.),  i  =  1,2,. . .  ,N 

~x.  ~x  1 

i 

for  cubic  approximations.  Here 


(2.25) 


Un  :=  U(xn,t  );  Un 
-l  ~  l  n  x. 


:=  U  (X.  ,t  ) 
x  i  n 


We  obtain  the  approximate  boundary  conditions  (2.6)  by  substi¬ 
tuting  (2.3)  into  (1.3),  integrating  the  resulting  equation  from  t  to  t^+1 ,  and 
evaluating  the  integrals  by  the  Trapezoidal  rule.  Each  boundary  condition  is 
associated  with  a  particular  partial  differential  equation  in  the  vector  system. 
The  test  space  1/  is  modified  by  setting  the  test  functions  4^  and  4N  (for  linear 
approximations)  or  $2n-1  ^or  cubic  approximations)  equal  to  zero  for  those 
partial  differential  equations  associated  with  boundary  conditions.  This  has 
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'C 


the  effect  of  replacing  the  Galerkin  approximation  of  a  partial  differential 
equation  at  either  x  =  0  or  x  =  a  by  its  corresponding  approximate  boundary 

condition. 


The  system  (2.6),  (2.23)  is  a  nonlinear  algebraic  system  for  determining 

n+1  ,  .  ,  ,  .  .  n+1  „n+l  ,  .  .. 

,  l  =  1,2,...,N  for  linear  approximations  or  U .  ,  U  ,  l  =  1,2,..., I.,  for 

i 

cubic  approximations  given  the  same  information  at  t  =  t^ .  We  solve  this  non¬ 
linear  system  by  Newton's  method  which  requires  the  computation  of  the  Jacobian 

of  the  vector  [FfU,^)  (FfU,^)  ,  —  ,F(U,4>K)]T  with  respect  to  [u^  1  ,V™  \...,U^  1]T 

,  .  .  .  .n+1  ,,n+l  „n+l  T,n+1  ,,n+l  n+l.T  c  w.  • 

for  linear  approximations  or  [U  ,  U  ,  U  ,  U  ,...,U  ,  U  ]  for  cubic 

-1  ~Xx  -X2  -N  -X15 

approximations.  The  Jacobian  will  be  block  tridiagonal  because  of  the  local 
nature  of  4..  The  elements  in  the  i  th  block  of  rows  will  be  the  M  x  M  matrices 

l 


(2.26a) 


3F(U,4»  ) 


3un+1 


j  =  i  -  1 ,  i,  i+1 


for  linear  approximations  and  the  2M  x  2M  matrices 


(2.26b) 


9r^2i-i) 

3un+1 

~3 


3F(u,4>2i) 

3un+1 

3 


3f  (u,<t>2i_i> 


3u 


n+1 

x . 

3 


,  j  =  i  -  1,  i,  i  +  1, 


3F(U,<f>2i> 


9u 


,n+l 


for  cubic  approximations.  The  elements  of  (2.16)  are  obtained  from  (2.22,23)  in 
a  relatively  straightforward  manner,  but  their  computation  requires  users  of  our 
code  to  supply  subroutines  that  define  f  (x,t,u,u  ),  f  (x,t,u,u  ),  and  D  (x,t,u). 

~u  —  -X  -U  -  n. x  -u  ~ 

-X 

Subroutines  that  define  f(x,t,u,u  )  and  D(x,t,u)  must,  of  course,  also  be  supplied. 

~  -  ~x  ~  - 

We  calculate  and  factor  the  Jacobian  once  per  time  step  and  use  Utx.t^)  as  an  initial 
guess  for  U(x,t  ,).  The  linearized  Newton  system  is  solved  by  an  efficient  block 
tridiagonal  algorithm  that  uses  pivoting  both  within  and  outside  of  blocks  (cf. 

Davis  [11]).  Generally  two  iterations  are  performed  per  time  step. 
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In  Section  2  we  developed  a  finite  element  method  to  obtain  numerical 


solutions  to  systems  of  partial  differential  equations  on  nonuniform  trape¬ 
zoidal  grids.  In  this  section  we  construct  an  algorithm  to  select  a  grid  at  . 

t  =  t  ,  so  that  the  L_  norm  of  the  error  at  t  ,  is  approximately  minimized. 
n+1  2  n+1 

This  algorithm  builds  upon  the  work  of  deBoor  [12]  ,  Lawson  [26] ,  and  Jupp 

[23]  on  variable  knot  spline  interpolation. 

For  most  of  this  section  we  will  be  discussing  approximations  at  a  single 

time  level,  say  t  =  fc  ,  so  whenever  there  is  no  possibility  of  confusion  we 

n 

omit  the  n  superscript  on  x?  and  U?  and  surpress  the  t  dependence  when  writing 
u(x,t).  We  also  present  the  development  for  scalar  functions  u  and  indicate 
the  extensions  to  vector  functions  in  Section  3.2. 

It  is  well  known  (cf.  [10,35,36])  that  the  errors  in  finite  element- 
Galerkin  methods  for  problems  like  (1.1-3)  satisfy  estimates  of  the  form 


(3.1)  l l u— U l (  <  cl lu-Pul I  , 

2  ~  ~  L2 

where  PU £  U  interpolates  u.  Thus,  the  error  in  the  solution  of  the  partial 

•»  Jv 

differential  equation  is  bounded  by  an  interpolation  error.  The  following 
result  (cf.  e.g.  Pereyra  and  Sewell  [31])  indicates  how  to  minimize  this 
interpolation  error  for  piecewise  polynomial  interpolants . 


Lemma:  Let  II  :=  (0  =  x.  <  x  <•••<  x  =  a)  be  a  partition  of  [0,a] 

N  12  N 

£+1 

into  N-l  subintervals  and  let  u(x)€C  [0,a]  .  The  piecewise  poly¬ 


nomial  of  degree  Z  on  i  =  If  2, - ,N-1,  that  interpolates 

to  u  on  nN  has  minimal  L^  error  when  the  knots  x^  ,  i  =  2, 3,..., N-l 
are  chosen  such  that 


(3.2)  h*+1|uU+1) (C±) I  =  E,  i  =  1,2 . N-l, 

(£) 

where  u  is  the  Z  the  derivative  of  u  with  respect  to  x,  e  (x.,x^+^),E 


is  a  constant,  and 
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(3.3) 


h .  = 
1 


'i+1 


x. . 

i 


The  Lemma  states  that  the  interpolation  error  is  minimized  by  selecting  the 
partition  in  such  a  way  that  the  quantity  h^+1 | u (£  ) |  is  equidistributed . 
Considerable  success  has  been  achieved  by  using  this  result  to  implement  adaptive 
grid  algorithms  for  two-point  boundary  value  problems  (cf.  Lentini  and  Pereyra 
[27],  Ascher,  Christiansen,  and  Russell  [2],  or  Russell  and  Christiansen  [34]). 
Nevertheless,  some  practical  difficulties  still  remain  and  we  discuss  these  and 
our  solutions  to  them  below. 

Rather  than  work  with  (3.2)  directly,  we  follow  Lawson  [27]  and  Jupp  [23] 
and  express  (3.2)  in  the  form 


(3.4)  p.  :=  h.+1/h.  =  [|u(Ul)(Ci)/ua+1)(C.+1)!]1/(Ul),  i  =  1,2 . N-2, 

where  2=1  for  piecewise  linear  and  2.  =  3  for  piecewise  cubic  approximations. 
In  addition  to  (3.4)  we  impose  the  constraint  that 


(3.5) 


hl  +  h2  + 


+  h. 


N-l 


V 


This  can  be  expressed  in  terms  of  the  p^'s  by  defining 

(3.6)  z  :=  1  +  (p^  +  (PXP2)  +  (P1P2P3)  +  •••  +  (P1P2P3  •••  Pn-25 
and  observing  that 

(3.7)  z  =  (hp  +  h2  +  h3  +  . . .  +  hN_1)/h1=  (*N  _  Xj.J/hj. 

Equations  (3. 3, 4, 6, 7)  permit  us  to  determine  h^,  i  =  1,2,..., N-l,  and 
x^,  i  =  1,2,...,N,  in  terms  of  u^+^  without  an  explicit  determination  of  E. 

Of  course,  u^+^  is  unknown  and  must  be  approximated  by  The 

finite  element  procedure  provides  us  with  an  approximate  solution  U  and  for 
cubics  an  approximate  first  derivative  Ux>  However,  equation  (3.4)  requires 
a  knowledge  of  second  derivatives  for  linear  approximations  and  fourth  deriva¬ 
tives  for  cubic  approximations.  This  situation  typically  arises  in  adaptive 
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mesh  algorithms  and  it  is  usually  resolved  by  using  finite  difference  approxima¬ 
tions  for  the  necessary  higher  derivatives. 


DeBoor  (12]  used  finite  difference  approximations  to  choose  mesh  points  for 
the  solution  of  two-point  boundary  value  problems  by  assuming  that  the  (P.-»l)st 
derivative  was  constant  on  each  subinterval.  We  modify  this  scheme  slightly  by 
assuming  that  is  linear  on  each  subinterval  and  takes  on  the  following 


values  at  the  nodes: 


(3.8a) 


where 


’4UV2/IW' 

2AUl/2/(W ' 


i  =  1 


l  = 


ItU+l)  ,  x  V 
U  (xJ 


iUi-3/2/fhi-l+hi-2>  ♦  «£U/<VW.  1  "  3’4 . 

4oiS-U/,Vi*Vj''  i‘K 


(3.8b) 

We  use  the  approximation 
(3.9a)  U 

for  linear  polynomials  and 


i-l/2  *  <Ui  - 


(3.9b) 


“i-l/2  "  1JlUl-rUi)/hi-l 


+  6 (U  +0  )/h.  , 

xi-l  *i  1‘1 


for  cubic  polynomials. 

We  note  that  p^  becomes  infinite  or  indeterminate  when  u^  ^  (t^)  =  0 
(cf.  (3.4));  hence*  we  can  expect  numerical  difficulties  when  u^  (x)  is  small 
on  any  subinterval.  Indeed  numerical  experiments  have  shown  that  the  mesh 
becomes  very  sensitive  to  small  perturbations  whenever  ^  (x)  is  of  the  same 
order  of  magnitude  as  the  discretization  error  in  the  computed  solution  U. 
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We  combat  this  problem  by  imposing  a  lower  bound  on  +  (x)  j.  Thus, 

we  let  At  *  (t  , -t  )  and  h  ■  a/N  denote  the  current  time  step  and  the  average 

n  n+x  n  ’ 

mesh  spacing,  respectively,  and  for  linear  approximations  we  calculate  'u  (x  )! 

i 

2  2  2 

as  the  maximum  of  the  value  computed  by  (3.8,9)  and  max  (At  /h  ,h  )  while  for 
cubic  approximations  we  calculate  |l/iv^(x^)|  as  the  maximum  of  the  value  com¬ 
puted  by  (3.8,9)  and  12  max(At/h,h2)  +  6  max  (At4/h3 ,h2) .  These  limits  were 
determined  empirically.  They  are  small  enough  so  that  they  do  not  affect  the 
mesh  adaption  procedure  when  (x)  is  not  small  but  large  enough  to  avoid  the 

numerical  difficulties  caused  by  vanishing  values  of  U^  +  1'(x).  Observe  that  if 
^  (x)  is  uniformly  small  on  [0,a]  our  limits  assure  that  the  solution  of 
(3. 3, 4, 6, 7)  is  a  uniform  mesh,  as  it  should  be  in  this  case. 

The  discussion  thus  far  has  concerned  the  computation  of  an  optimal  grid 
at  a  time  level  t  where  the  solution  un  has  already  been  computed.  Ke  wculd 
also  like  to  estimate  an  optimal  grid  at  time  level  t^+^  prior  to  computing  the 
solution  there.  This  can  be  done  by  extrapolating  the  optimal  grids  computed 
at  a  number  of  previous  time  levels  to  t  +^.  It  was  somewhat  surprising  that 
numerical  experiments  seemed  to  favor  zero  order  extrapolation;  i.e.,  the  optimal 

grid  computed  at  time  level  t  is  used  at  time  level  t  . .  Multi-level  extrapo- 

n  n+1 

lation  consistently  overestimated  the  distance  that  a  mesh  point  should  move  in 
one  time  step  and  then  overcorrected  this  error  in  the  next  time  step.  Ir.  some 
cases  this  caused  the  mesh  to  oscillate  wildly  when  in  fact  the  solution  changed 
very  little.  When  we  simply  extrapolated  the  optimal  mesh  determined  at  the  prev¬ 
ious  time  level  it  tended  to  follow  the  solution  even  when  rapidly  moving  fronts 
were  present. 

It  is  easy  to  Bhow  that  the  mesh  selection  strategy  (3. 3, 4,5,b)  maintains 
the  knot  ordering  so  that  no  two  mesh  points  can  cross.  It  does  not  however 
prohibit  severely  distorted  trapezoidal  elements.  Ciarlet  and  Raviart  [10]  and 
Babuska  and  Aziz  [3]  have  studied  the  effect  of  element  distortion  on  the 
accuracy  of  the  finite  element  method.  They  have  shown  that  the  error  obtained 
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when  computing  on  trapezoidal  elements  is  a  multiple  of  the  error  obtained  when 
computing  on  rectangular  elements.  The  multiplicative  factor  is  proportional  to 
a  power  of  the  magnitude  of  the  derivatives  of  the  transformation  (2.8) .  There¬ 
fore  we  must  control  the  magnitude  of  these  derivatives  in  order  to  maintain 
acceptable  accuracy.  We  let 

(3.10)  h?  -  x"  -  x"  At  =  t  .  -  t  .  tan  u.  -  hn/At  . 

i  i+l  i  n  n+1  n  1  1  n 

Hence,  ^  is  the  angle  between  the  line  x^  (t)  and  the  positive  t  axis. 

Differentiating  (2.8)  and  using  (3.10)  we  find 
-  [hj  +  T(h"+1-h^)]/2, 

(3.11)  Xt  *  At^[tan  ok  +  (tan  <oi  +  1-tan  <o^)  (£+l)/2], 

t_  *  0,  t  *  At  . 

5  T  n 

Since  the  magnitudes  of  h?  and  h”+^  are  controlled  by  the  bounds  that  we 
imposed  on  |u^+^|  and  At^  is  prescribed,  we  can  limit  the  magnitude  of  the 
derivatives  in  (3.11)  by  controlling  the  growth  of  |tan  on  | .  We  found  that  the 
condition 

(3.12)  max  (oj  |  <_  3tf/8 

l<i<N  i 

worked  well  in  practice. 

3.1.  Mesh  Selection  Algorithm 

In  this  section  we  discuss  some  details  of  a  mesh  selection  algorithm  based 

on  the  discussion  of  the  previous  section.  The  first  algorithm  uses  the  finite 

element  solution  (!(x,t  )  calculated  on  the  mesh  x?,  i  ■  1,2,...,N,  and  equations 

n  1 

(3. 3, 4, 6, 7)  to  find  a  new  mesh  at  t  ■  t  that  satisfies  the  optimality  condition 

n 

(3.2).  This  is  the  mesh  that  should  have  been  used  to  calculate  U(x,t  ).  Instead 

n 

we  use  it  in  the  second  algorithm  to  estimate  an  optimal  mesh  at  t  =  t  . 
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The  difficulty  in  solving  (3. 3, 4, 6, 7)  for  the  optimal  mesh  is  that  these 
equations  are  nonlinear  and  must  be  solved  iteratively.  We  use  a  relaxation 
scheme  that  is  similar  to  one  which  has  been  analyzed  by  Isaacson  and  Keller 
[21,  Chap.  3].  They  give  necessary  convergence  criteria,  but,  we  chose  not 
to  incorporate  these  into  our  algorithm  because  they  require  too  much  additional 
computation.  The  following  algorithm,  which  calculates  the  relaxation  parameter 
heuristically,  has  not  failed  to  converge  in  any  of  our  tests. 

1.  Set  the  relaxation  parameter  Q  -.=  1  and  let  x|^  ••  = 

x?,  i  =  1,2,...,N  be  an  initial  guess  for  the  optimal 
mesh.  Calculate 


,(0)  :=  <X<0)  -X^J/Cxf-X^) 


(1) 


:=  z  +  2€ 


V  :  =  1 

where  €  is  a  convergence  tolerance. 

2 .  Compute 

„a+l)  (0),  .  .  , 

U  (x^  ) ,  i  -  1,2 , . , . ,N 

using  equations  (3.8,9). 

3.  While  | z  -z  ^V+1 M  >  €  or  V  <  V  do 

-  '  —  max  — 

4.  Calculate  (x|V  ^),  i  =  1,2,...,n  by  linear 

interpolation  of  U^+1^(x|°^),  i  =  1,2,...,N. 
Calculate 

i  *  1,2, .. . ,N-2. 


and 

:  (v) 


,  ,  (V),  ,  (V)  (V). 

1  +  (Pi  )  ♦  (Pj  P2  )  + 


,  (V)  (V)  (v) 

•  (P1  P2  'PN-2 
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5.  If  v  >  1  then 


if  .;<v)_2cv-x,| >  li(v-i).8(v-2)Ithenn  :.n/2 


6.  Calculate 


(V)  ,  (v-l) 


1 

(v) 


(xn'  "  *iV’1))/z(V> 


(V)  (V-l) 

X1  :  =  X1 


(V) 


(V) 


(v-l) 


(V) 

._  ~<v> 

+ 

i+l 

■ —  X  . 

1 

1 

(v) 

‘i+l 

h.(V) 

l 

,(v) 

‘i+l 

~(v) 

:=ftci+i 

* 

,(V) 

.  ,  (V) 

-  (XN 

-  x^v))/(x 

i  =  1,2, ,  N-2 


7.  V  :*  v  +  1 

(v) 

For  vector  systems  we  need  only  to  change  the  definition  of  p^  used  in  step  3. 


We  used 


where  IP  is  the  j  th  component  of  U. 

After  we  compute  a  convergent  mesh,  x^+1  =  x|V^  '  * 


1,2,...,K,  we  perform 


the  following: 

1 .  Compute 


.  t*n+l  ni 

Ax  “  max  x .  -x . 
max  2<i<N-l  *  1 


2.  If  Ax  <  At  tan  (3it/8) 
—  max  —  n 


then  Ax,,  !*  Ax 
-  fix  max 


else  Ax,,  !«  At_  tan(3it/8). 
fix  ft 
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H 


Compute  corrected  mesh  xn+^  as 


x?+1  :=  xn  +  (x?+1-x?)Ax  /Ax  .  , 
1  111  max  fix 


i  =  2,3,... , N-l 


Steps  2-3  prevent  the  elements  from  becoming  too  distorted. 

The  algorithms  contain  several  approximations  and  heuristic  procedures. 

Derivatives  are  estimated  by  differences  and  are  assumed  to  vary  linearly 

between  mesh  points.  Zero  order  extrapolation  was  used  to  predict  optimal 

grids  at  subsequent  time  levels.  Grids  were  restrained  to  prevent  severe 

element  distortion.  Even  with  these  approximations  the  mesh  selection 

algorithms  performed  satisfactorily  on  all  test  examples  that  we  considered. 

In  addition  we  note  that  Rheinboldt  [32]  has  shown  that  an  order  A  error  in 

2 

the  placement  of  the  optimal  mesh  only  produces  an  order  A  change  in  the 
computed  solution.  Thus,  it  suffices  to  only  be  close  to  the  optimal  mesh 
in  order  to  reap  its  benefits. 


-18- 


4.  Computational  Results 

In  this  section  we  examine  the  performance  of  our  method  on  four  problems 
which  are  graded  in  difficulty  such  that  each  one  exercises  an  additional  facet 
of  the  method.  The  following  norms  are  used  to  evaluate  the  performance  of  our 
method  on  examples  where  exact  solutions  are  known. 

(4.1a)  I  I e  (t) | |  :=  max  |e(x.,t)|  :=  max  |u(x. ,t)  -  U(x. ,t)| 

'  00  1  i  'CO  1  -  1  00 

l<i<N  l<i<N 

2  N-i  , 

(4.1b)  ||e(t)|iL  :=  l  (h./2)(|e(xi,t)|oo  +  |e(x.+1>t)|j, 

2  i=l 

where 


(4.1c) 


M.  :=  max  |v 
l<k<M 


Example  1: 


“t  ■  (1/”  "xx-  0  <  »  <  1.  t  >  o 


(4.3) 


u (x ,0)  =  sin  tt  x,  u(0,t)  =  u(l,t)  =  0. 
The  exact  solution  is 


u(xft)  =  e  t  sin  it  x. 

Analysis  presented  in  [11]  indicates  that  the  finite  element  method  des- 

2  2 

cribed  in  section  2  would  have  error  of  0(h)  +  0(At  )  with  linear  elements 
4  2 

and  0(h)  +  0(At  )  with  cubic  elements  on  a  uniform  spacial  mesh  of  width  h  and 
a  uniform  time  step  of  duration  At.  We  created  this  simple  constant  coefficient 
example  to  verify  that  these  errors  are  actually  attained.  Figures  2  and  3 
present  plots  of  the  error  at  t  =  1  as  a  function  of  h  for  linear  and  cubic 
approximations,  respectively. 

The  analysis  of  [11]  predicts  that  the  points  on  Figure  2  for  which  At  =  h 

2 

and  the  points  in  Figure  3  for  which  At  *  h  should  lie  on  straight  lines  having 
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I 


slopes  2  and  4,  respectively.  These  lines  are  shown  confirming  that  the  theor¬ 
etical  error  bounds  are  actually  attained. 

Example  2 : 

(4.3a)  u  =  Ou  +  f (x) ,  0  <  x  <  1,  t  >  0,  O  >  0. 

t  XX 

The  initial  conditions,  boundary  conditions,  and  source  f  are  chosen  so  that  the 
exact  solution  is 

(4.3b)  u(x,t)  =  tanh  (r^  (x-1)  +  r  t) 

The  solution  (4.3b)  is  a  wave  that  travels  in  the  negative  x  direction 
when  r^  and  r^  are  positive.  The  values  r^  and  r^  determine  the  steepness  of 
the  wave  and  its  speed  of  propagation.  Thus,  the  problem  can  be  made  more  or 
less  difficult  by  adjusting  r^  and  r^. 

We  created  this  problem  to  study  the  effectiveness  of  our  adaptive  mesh 
algorithm  at  concentrating  grid  points  in  transition  regions,  following  moving 
fronts,  and  reducing  errors  below  those  of  uniform  grid  calculations. 

We  first  solve  problem  (4.2)  with  r^  =  r^  =  5,  uniform  time  steps  of 
At  =  0.01,  10  elements  per  time  step,  and  linear  approximations.  The  mesh  com¬ 
puted  by  our  adaptive  mesh  algorithm  is  shown  in  Figure  4.  The  grid  points  are 
concentrated  in  the  region  of  maximum  curvature  and  move  to  the  left  with  the 

wave.  As  the  wave  front  passes  out  of  the  domain  and  u  becomes  small,  the 

xx 

grid  points  move  toward  a  uniform  distribution.  It  is  clear  that  the  grid 
adapts  to  the  solution  and  follows  its  progress. 

As  a  somewhat  more  difficult  problem  we  solve  (4.3a)  with  initial  con¬ 
ditions,  boundary  conditions,  and  forcing  function  chosen  so  that  the  solution 
is  given  by  (4.3b)  with  r^  *  r^  =  100.  The  wave  front  is  much  steeper  than  in 
the  previous  test  of  (4.3). 
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In  Table  1  we  present  a  comparison  at  t  =  1  of  the  results  of  computation 
using  linear  approximations  on  a  variety  of  uniform  and  variably  spaced  meshes. 
These  results  are  somewhat  disappointing.  At  best  the  mesh  moving  scheme  improves 
the  accuracy  of  the  solution  only  slightly.  The  improvement  is  greatest  when  it 
is  small  and  in  some  cases,  when  At  is  large,  the  uniform  mesh  is  more  accurate. 

A  closer  examination  explains  these  results  and  reveals  something  about  the  nature 
of  this  mesh  moving  scheme. 

Table  1  shows  that  the  solution  of  this  problem  was  not  computed  accur¬ 
ately  with  either  a  uniform  or  a  variable  mesh.  This  can  be  explained  by  exar.in- 

* 

ing  the  time  evolution  of  the  solution  at  a  fixed  value  of  x,  say  x  .  The 

solution  is  approximately  given  by  -1  until  the  time  when  the  wave  reaches  the 

* 

point  x  .  It  then  jumps  suddenly  to  a  value  near  1.  If  the  time  step  At  is 
too  large  to  resolve  this  transition,  we  would  expect  large  errors  in  the 
vicinity  of  the  wave.  The  solid  curve  in  Figure  5  confirms  this  prediction. 

The  mesh  selection  procedure  misinterprets  these  errors  as  being  part  of  the 
solution  and  places  too  many  points  in  the  region  outside  of  the  wave  front. 

Thus,  a  suboptimal  mesh  is  selected  and  the  expected  decrease  in  the  error  is 
not  obtained.  When  At  becomes  small  enough  to  adequately  resolve  the  passing 
wave  the  mesh  selection  procedure  does  improve  the  accuracy  of  the  solution 
(cf .  Table  1) . 

This  points  out  the  need  for  an  algorithm  to  adaptively  refine  time  steps 
in  the  vicinity  of  severe  temporal  gradients.  Such  a  procedure  was  used  by 
Berger  et  al.  [5]  to  solve  hyperbolic  partial  differential  equations  and  we  are 
currently  studying  its  suitability  for  our  code. 

Table  2  summarizes  the  results  of  computations  performed  on  the  same 
problem  using  cubic  approximations.  In  these  cases  the  time  steps  At  were 
small  enough  to  resolve  the  transition  of  the  solution  and  the  cubic 


-21- 


approximations  were  accurate  enough  to  provide  us  with  reasonable  estimates  of 
the  derivatives.  As  a  result,  the  variable  mesh  scheme  improved  the  solution 
significantly. 

Figures  5  and  6  for  linear  and  cubic  approximations,  respectively,  show 
that  the  mesh  selection  algorithm  tends  to  distribute  the  local  error  evenly 
over  the  domain  and  thus,  as  indicated  in  Section  3,  approximately  minimizes 
the  error  in  L^. 

Example  3:  (Burgers'  Equation) 


(4.4) 


u  =  -uu  +  Cu  ,  0  <  x  <  1,  t>0, 

t  X  XX 


and  C 


5  x  10“3. 


u(x,0)  =  sin  7T  x,  u(0,t)  =  u(l,t)  =  0, 


It  is  well  known  that  the  solution  to  this  problem  is  a  wave  that  steepens 
and  moves  to  the  right  until  a  shock  layer  forms  at  x  =  1 .  After  a  time  of  0(1/0 
the  wave  dissipates  and  the  solution  decays  to  zero.  Figures  7  and  8  show  the 
results  of  computations  on  this  problem  using  linear  approximations  on  a  uniform 
mesh  and  a  variable  mesh  with  a  constant  time  step  of  At  =  0.1  and  10  elements 
per  time  step.  The  results  in  Figure  7  are  typical  of  finite  difference  or  finite 
element  calculations  for  this  problem.  Spurious  oscillations  develop  in  the  com¬ 
puted  solution  unless  the  mesh  width  is  of  the  same  order  as  the  width  of  the 
shock  layer,  which  is  0(/F)  for  this  example.  The  variable  mesh  results  in  Figure 
8  largely  surpress  these  oscillations  by  automatically  concentrating  the  mesh  in 
the  shock  region  as  the  wave  steepens. 

When  Example  3  is  solved  using  cubic  approximations  on  a  uniform  mesh  we 

find  that  the  solution  U?  at  the  nodes  is  computed  accurately;  however,  there  are 

large  errors  in  the  slope  of  the  solution  U°  at  the  nodes  when  the  mesh  is  not 

xi 

suitably  fine  in  the  shock  region.  This  effect  is  exhibited  in  Figure  9  where 
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the  solution  at  t  =  0.6  is  shown  for  a  calculation  performed  with  At  =  0.1  and 
N  =  10.  Equations  (2.6,19,20,23)  were  used  to  calculate  the  solution  between 
mesh  points. 

One  possible  explanation  of  this  behavior  was  proposed  by  Killer  and 
Miller  [29] ,  but  they  do  not  explain  why  the  large  error  in  the  slopes  do  not 
feed  back  and  cause  large  errors  in  the  function  values. 

Once  again,  these  problems  are  corrected  when  the  mesh  adapts  with  the 
solution.  Figure  10  shows  the  result  of  a  similar  computation  using  cubic 
approximations  on  a  variable  mesh.  Both  the  function  values  and  slope  values 
are  computed  accurately  at  the  nodes. 

Example  4: 


(4.5) 


b  =  [y  (s)b  ]  -  [by  (s)  s  ]  , 

t  XX  X  X 


s  =  -k  (s) b ,  0  <  X  <  5,  t  >  0. 

This  two  component  nonlinear  system  was  studied  by  Keller  and  Odell  [24,30]  as 

a  model  for  the  chemotactic  motion  of  bacteria.  The  quantity  b(x,t)  denotes 

the  bacterial  density  and  s(x,t)  denotes  the  concentration  of  the  critical 

substrate  (bacterial  food).  If  the  functions  y,X  and  k  satisfy  conditions 

derived  by  Keller  and  Odell  [25],  equations  (4.5)  have  travelling  wave  solutions. 

These  solutions  have  been  computed  by  Odell  and  Keller  [30]  and  are  interpreted 

as  travelling  bands  of  bacteria.  For  our  study  we  choose  k(s)  =  1,  y (s)  =  y  , 

and  Y(s)  =  6  /s  where  y  and  6  are  constants.  The  initial  conditions  are 
A  o  o  o 

shown  in  Figure  11a  and  the  boundary  conditions  are 


(4.6) 


b (0,t)  =  b (5 ,t)  =  0,  s  (0,t)  =  1. 


We  solved  this  problem  for  6  /y  =2  using  cubic  approximations,  uniform 

o  o 

time  steps  of  At  =  0.005,  and  50  elements  per  time  step.  The  computed  solutions 
at  t  ■  0,  0.1,  0.5,  and  1.0  are  shown  in  Figures  11a, b,c,  and  d,  respectively. 
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The  method  places  the  majority  of  the  mesh  points  in  the  regions  of  the  wave 


fronts  and  follows  the  bacterial  motion.  The  results  indicate  that  our  adap¬ 
tive  mesh  algorithm  may  be  also  used  for  vector  systems  of  equations. 
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5.  Discussion  and  Conclusions 

The  computations  presented  in  the  last  section  show  that  it  is  possible  to 
construct  an  accurate  and  stable  adaptive  grid  finite  element  method  for  nonlinear 
systems  of  partial  differential  equations  and  that  such  techniques  offer  advantages 
ver  fixed  grid  techniques.  In  particular,  we  have  shown  that  the  error  estimates 
obtained  by  Davis  [11]  are  actually  realized  in  practice  and  that  the  adaptive  mesh 
algorithm  correctly  concentrates  the  mesh  in  a  sharp  transition  and  is  able  to 
follow  moving  fronts.  Examples  3  and  4  of  Section  4  indicate  that  our  method  is 
also  useful  for  nonlinear  equations  and  vector  systems  of  equations. 

In  the  present  study  we  used  piecewise  polynomial  functions  for  both  the 
trial  and  test  spaces.  However,  recent  work  of  Flaherty  and  Mathon  [15]  ,  Heinrich 
et  al.  [18] ,  and  Hemker  [19]  indicates  that  exponential  and  "upwinded"  polynomial 
functions  may  give  superior  test  spaces  for  singularly  perturbed  problems.  We 
plan  to  incorporate  these  functions  into  our  methods  shortly. 

All  of  our  calculations  were  performed  with  a  constant  time  step.  Examples 
3  and  4  of  Section  4  indicate  that  it  would  be  most  desirable  to  be  able  to  vary 
the  time  step  during  the  calculation.  Our  code  presently  allows  for  this,  but  as 
yet  we  have  not  implemented  an  algorithm  to  adaptively  alter  the  time  step.  We 
also  plan  to  add  this  feature  to  our  code  shortly. 

Other  areas  for  future  study  include  free  boundary  problems  and  higher 
dimensional  problems.  The  present  work  has  shown  that  it  is  possible  to  construct 
a  practical  adaptive  grid  finite  element  method.  Future  work  must  refine  this 
method  and  apply  it  and  test  it  on  a  greater  variety  of  problems. 
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TABLE  1 


Results  of  Computations  at  t  *  1  for  Example  2  with  r^  *  r  =  i0°  usin9  Linear 
Approximations  on  Uniform  and  Variably  Spaced  Grids. 


N 

At 

Uniform  Spacing 

Variable  Spacing 

1 lei It 

L2 

mm 

11*11. 

L2 

IUIL 

10 

0.168 

0.137 

0.459 

1.346 

1.107 

1.708 

0.492 

0.949 

0.146 

0.254 

0.121 

0.340 

20 

0.365 

1.391 

0.567 

1.00 

0.05 

0.177 

0.392 

0.155 

0.746 

0.01 

0.768  E-l 

0.226 

0.166  E-l 

0.870  E-l 

40 

0.025 

0.367  E-l 

0.697  E-l 

0.348  E-l 

0.565  E-l 

0.01 

0.342  E-l 

0.158 

0.106  E-l 

0.105 

100 

0.01 

0.701  E-2 

0.703  E-2 

0.493  E-2 

0.158  E-l 

0.144  E-2 

0.275  E-2 

TABLE  2 


O  At  =  0.05 
□  = 0.025 

A  =0.01 
0  =  0.005 


1.0  x  10 


At  =  0.5 


1.0  x  10' 


1.0  x  10  f| 


.05  .10 


.50  1.0 


Figure  3j  error  vs.  h  for  Example  1  computed  on  uniform  meshes  with  cubic 
approximations.  The  dotted  line  connects  points  for  which  At  »  h2. 


Figure  5s  Local  error  at  t  «  1.0  for  Example  2  with  r^  *  r2  =  100,  uniform  time 
steps  of  Lt  ■  0.01,  N  ■  20,  and  linear  approximations.  The  solid  curve  was  com¬ 
puted  on  a  fixed  uniform  mesh,  the  broken  curve  on  a  variable  mesh. 


Figure  6:  Local  error  at  t  =  1.0  for  Example  2  with  r^  «=  =  100,  uniform  time 
steps  of  At  *=  0.0025,  N  ■  20,  and  cubic  approximations.  The  solid  curve  was  com 
puted  on  a  fixed  uniform  mesh,  the  broken  curve  on  a  variable  mesh. 


Figure  6s  Local  error  at  t  «=  1.0  for  Example  2  with  r^  *  r^  *=  100,  uniform  time 
steps  of  At  =  0.0025,  N  >=  20,  and  cubic  approximations.  The  solid  curve  was  com¬ 
puted  on  a  fixed  uniform  mesh,  the  broken  curve  on  a  variable  mesh. 
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Lgure  7:  Solut 


ions 


ion  of  Example  3  for  various  values  of  t  using  linear  approxima- 
rm  mesh  with  At  •  0.1  and  N  ■  10. 


\  of  Example  3  for 
of  At  ■  0.01,  and 


Figure  11:  Computational  results  for  Example  4  with  ^Q/VQ  “2.0  using  uniform  time 
steps  of  At  “  0.005,  N  -  50,  and  cubic  approximations  at  t  ■  0,  0.1,  0.5,  and  1.0. 


